library(Seurat)
library(tidyverse)
library(tidyseurat)
library(EnhancedVolcano)
library(ggpubr)

# MS-CSF ------
sobj <- read_rds('DE_cells/data/ms_sobj.rds')
sobj |> colnames() |> head()

zxk_in_ms <- read_csv('DE_cells/results/barcodes_MS-CSF-DE.csv') |>
  set_names(c('.cell', 'DE.cells'))
zxk_cd3_ms <- read_csv('DE_cells/results/barcodes_MS-CSF-CD3-DE.csv') |>
  set_names(c('.cell', 'CD3neg.DE.cells'))

sobj <- sobj |>
  left_join(zxk_in_ms) |>
  left_join(zxk_cd3_ms)

csf_memory <- sobj |>
  filter(!is.na(DE.cells) & str_detect(monaco_label, 'Switched'))

zxk_degs <- csf_memory |>
  FindMarkers(ident.1 = 'TRUE', group.by = 'DE.cells') |>
  rownames_to_column('gene') |>
  as_tibble()

zxk_degs |>
  ggplot(aes(avg_log2FC, -log(p_val_adj))) +
  geom_point()

EnhancedVolcano(zxk_degs,
                lab = zxk_degs$gene,
                x = 'avg_log2FC',
                y = 'p_val_adj',
                title = 'Differential expressed genes in DE switched memory B cells from CSF',
                titleLabSize = 15,
                subtitleLabSize = 0,
                pCutoff = .001,
                drawConnectors = TRUE)

library(enrichR)
enrichR::listEnrichrDbs() |>
  filter(str_detect(libraryName, 'Biological'))

Idents(csf_memory) <- 'DE.cells'

DEenrichRPlot(csf_memory,ident.1 = 'TRUE',enrich.database = 'GO_Biological_Process_2023', max.genes = 2000)

### manual enrichr
man_enrichr <- filter(zxk_degs, avg_log2FC > 1, p_val_adj < .05) |>
  pull(gene) |>
  enrichr('GO_Biological_Process_2023') |>
  pluck(1) |>
  as_tibble() |>
  filter(Adjusted.P.value < .05)

man_enrichr |>
  slice_head(n = 10) |>
  mutate(Term = fct_reorder(Term, Combined.Score)) |>
  ggplot(aes(Term, Combined.Score, fill = Adjusted.P.value)) +
  geom_col() +
  coord_flip() +
  theme_pubr() +
  scale_fill_viridis_c()

zxk_cd3neg.degs <- csf_memory |>
  FindMarkers(ident.1 = '1', group.by = 'CD3neg.DE.cells') |>
  rownames_to_column('gene') |>
  as_tibble()

zxk_cd3neg.degs |>
  ggplot(aes(avg_log2FC, -log(p_val_adj))) +
  geom_point()

csf_memory |>
  VlnPlot('nFeature_RNA', group.by = 'DE.cells') +
  stat_compare_means(method = 't.test')

csf_memory |>
  VlnPlot('nCount_RNA', group.by = 'CD3neg.DE.cells') +
  stat_compare_means(method = 't.test')
